
*************************************************
* Private Equity NH Ownership Time Trends - Figure 1
*************************************************
// Aggregate Graphs
use "$Data/Temp_Files/Aggregate_Base", replace
keep if year < 2020
replace capital_invested = capital_invested/1000 //make into billions
gen ind = 1 if year >= 2016
replace ind = 0 if missing(ind)

bgshade year, shaders(ind)  sstyle(noextend lcolor(gs14)) ///
twoway (bar capital_invested year if type=="all healthcare", yaxis(2) color("emidblue") || ///
	connected deal_count year if type=="all healthcare" , ///
	graphregion(color(white)) ///
	ytitle("Number of Deals", axis(1)) ytitle("Dollar Amount (Billions)", axis(2)) ///
	xtitle("Year") legend(order(1 "Capital Invested" 2 "Deal Count")) xlabel(2000(2)2020))
gr export "$Output/Figures/1A_All_Healthcare.png", replace
	
bgshade year, shaders(ind)  sstyle(noextend lcolor(gs14)) ///
twoway (bar capital_invested year if type=="elder and disabled care", yaxis(2) color("emidblue") || ///
	connected deal_count year if type=="elder and disabled care" , ///
	graphregion(color(white)) ///
	ytitle("Number of Deals", axis(1)) ytitle("Dollar Amount (Billions)", axis(2)) ///
	xtitle("Year") legend(order(1 "Capital Invested" 2 "Deal Count")) xlabel(2004(2)2020))
gr export "$Output/Figures/1B_Elder.png", replace


//Plotting Nursing Homes Deals
use "$Data/Temp_Files/Facility_Base", replace
// Keep relevant years
keep if year < = 2015
keep if year >= 2004
// Generate total number of patients
gen totpat = occupied_bed_pct * totbeds/100
gen all_firms = 1 // Indicator for all firms
gen pe_totpat = peind*totpat // PE patients
// Collapse
collapse (sum) all_firms pe_target = peind totpat pe_totpat,by(year) // Collapse by year
gen per_firm = pe_target / all_firms *100 
gen per_pat = pe_totpat / totpat * 100

//Deals 
twoway connected pe_target year,yaxis(1) ylabel(0(400)1600) ytick(0(200)1600) ytitle("# of Facilities",size(large)) lcolor("dknavy") || ///
 connected pe_totpat year,yaxis(2) ylabel(0(30000)150000,axis(2)) ytick(0(15000)150000,axis(2)) ytitle("# of Patients" ,axis(2) size(large)) lpattern("-") lcolor(" maroon") ///
 bgcolor(white) graphregion(color(white)) xtitle("Year",size(large))  xlabel(2004(2)2015) xtick(2004(1)2015) ///
 legend(label(1 "# of Facilities - PE Owned") label(2 "# of Patients - PE Owned") bmargin(tiny) keygap(1) symxsize(10)) 
// Save
graph export "$Output/Figures/1C_PE_deals_all.png", replace

// Percent Deals
twoway connected per_firm year, lcolor("dknavy") || ///
 connected per_pat year, lpattern("-") lcolor(" maroon") ///
 bgcolor(white) graphregion(color(white)) xtitle("Year" ,size(large))  xlabel(2004(2)2015) xtick(2004(1)2015) ///
 legend(label(1 "% of Facilities - PE Owned") label(2 "% of Patients - PE Owned") bmargin(tiny) keygap(1) symxsize(10)) ytitle("% PE - Owned" ,size(large)) ylabel(0(2)12) 
// Save
graph export "$Output/Figures/1D_PE_deals_pct_all.png", replace

*************************************************
* Initial Patient Assesments - Event Study - Figure 2
*************************************************
// Load Data
use "$Data/Temp_Files/Person_Base", replace

// Generate event time as Year relative to PE acquisition
destring deal_year, replace
gen event_time2 = 	year-deal_year
*** Event study dummies
replace event_time2 = 0 if mi(event_time2)
tab event_time2, gen(deal_dummies)
//buyout year is deal dummy 6

label var ADL_self_initial  "All ADL"
label var alzheimers_initial  "Alzhiemers"
label var dementia_initial "Dementia"
label var depression_initial "Depression"

eststo clear
estimates clear
* Output Vector
local mds_initial ADL_self_initial  alzheimers_initial dementia_initial depression_initial 
*Control vector;
global mktcntrl la_herfbeds_cty res_profit_cty res_multifac_cty res_hospbase_cty logcount agg_cmi_w agg_black_w acuityindx_w

global dummies deal_dummies1-deal_dummies4 deal_dummies6-deal_dummies11

foreach depvar of local mds_initial  {

	local vtext: variable label `depvar'

	* Regression
	reghdfe `depvar' $dummies $mktcntrl, absorb(year facid) tol(1e-6) accel(sd) cluster(facid)

	* Save for figure
	local var outcome
	gen `var'_b = .
	gen `var'_sd = .
	gen `var'_u = .
	gen `var'_l = .
	
	
	
	* pick years before; 
		forvalues x=2(1)4{
		local y = `x' - 1
		replace `var'_b=_b[deal_dummies`x'] in `y'
		replace `var'_sd=_se[deal_dummies`x'] in `y'
		replace `var'_u=_b[deal_dummies`x']+1.96*_se[deal_dummies`x'] in `y'
		replace `var'_l=_b[deal_dummies`x']-1.96*_se[deal_dummies`x'] in `y'
	}
	
	* Baseline
	replace `var'_b = 0 in 4 
	
	* fill in for years after 
	* UPDATE if data changes (i.e. deal year == 0 dummy can change)
	forvalues x=6(1)11{
		local y = `x' - 1
		replace `var'_b=_b[deal_dummies`x'] in `y'
		replace `var'_sd=_se[deal_dummies`x'] in `y'
		replace `var'_u=_b[deal_dummies`x']+1.96*_se[deal_dummies`x'] in `y'
		replace `var'_l=_b[deal_dummies`x']-1.96*_se[deal_dummies`x'] in `y'
	}
	
	* Figure label cleanup
	gen n=_n
	label define betas   1 "-4" 2 "-3" 3 "-2" 4 "-1" 5 "0" 6 "1" 7 "2" 8 "3" 9 "4"  , replace 
	label values n betas
		
	
	twoway (rcap `var'_u `var'_l n, color(gs11)) (scatter `var'_b n, msize(medium) mcolor(navy) msymbol(circle) )  if n <=9,  graphregion(color(white)) legend(off)  xlabel(1(1)9,valuelabel)  ytitle("`vtext'", size(medsmall))  xtitle("Years from Deal (0 is Deal Year)")   yline(0, lcolor(black)) xline(4.5, lcolor(red) lpattern(dash))
	graph export "$Output/Figures/2_`depvar'_fac_did.png", replace 
 
	drop `var'_b
	drop `var'_sd
	drop `var'_u
	drop `var'_l
	drop n 
}



*************************************************
* Differential Distance - Figure 3
*************************************************
// Binscatters
// Probability of going to PE and Differential Distance_Instrument
use "$Data/Temp_Files/Person_Base", replace

reghdfe peind c.age#i.sex dual, absorb(married race facid gp_old* bene_hrr#year) tol(1e-6) accel(sd) cluster(facid) residuals(pe_res)
reghdfe diff_dist_w c.age#i.sex  dual, absorb(married race facid gp_old* bene_hrr#year) tol(1e-6) accel(sd) cluster(facid) residuals(dist_res)
reg pe_res dist_res
local b = string(_b[dist_res]*10,"%9.3f")
local se = string(_se[dist_res]*10,"%9.3f")
binscatter pe_res dist_res, ytitle("Stay at PE Nursing Home" , size(large)) xtitle("Differential Distance to PE NH (Miles)" , size(large)) line(qfit) note("Slope (per 10 Miles): `b' (`se')") n(10) ///
savegraph("$Output/Figures/3_Binscatter_PE_FacFE.png") replace

// Charleson Score
reghdfe charles_old, absorb( c.age#i.sex married race facid bene_hrr#year dual) tol(1e-6) accel(sd) cluster(facid) residuals(charles_res)
drop dist_res
reghdfe diff_dist_w, absorb( c.age#i.sex married race facid bene_hrr#year dual) tol(1e-6) accel(sd) cluster(facid) residuals(dist_res)
reg charles_res dist_res
local b = string(_b[dist_res]*10,"%9.3f")
local se = string(_se[dist_res]*10,"%9.3f")
binscatter charles_res dist_res, ytitle("Differential Share High Risk Patients" , size(large)) xtitle("Differential Distance to PE NH (Miles)" , size(large)) note("Slope (per 10 Miles): `b' (`se')") yscale (range(-.05 .05)) n(10) ///
ylabel(-.05(.01).05) savegraph("$Output/Figures/3_Binscatter_charles_FacFE.png") replace


*************************************************
* Patient Outcomes - Differential Distance - Figure 4
*************************************************
/// Try to do binscatters
use "$Data/Temp_Files/Person_Base", replace
label var mortality_90  "Mortality (Stay + 30 Days)"
label var log_pmt  "Log Amount Billed per Patient Stay"
label var log_payment_total_post  "Log Amount Billed per Patient Stay + 90 Days"

reghdfe mortality_90 c.age#i.sex dual, absorb(married race facid gp_old* bene_hrr#year) tol(1e-6) accel(sd) cluster(facid) residuals(mort_res)
reghdfe log_pmt c.age#i.sex dual, absorb(married race facid gp_old* bene_hrr#year) tol(1e-6) accel(sd) cluster(facid) residuals(pay_res)
reghdfe log_payment_total_post c.age#i.sex dual, absorb(married race facid gp_old* bene_hrr#year) tol(1e-6) accel(sd) cluster(facid) residuals(pay_res_post)
reghdfe diff_dist_w c.age#i.sex  dual, absorb(married race facid gp_old* bene_hrr#year) tol(1e-6) accel(sd) cluster(facid) residuals(dist_res)

// BS 1
reg mort_res dist_res
local b = string(_b[dist_res]*10,"%9.3f")
local se = string(_se[dist_res]*10,"%9.3f")
binscatter mort_res dist_res, ytitle("Mortality (Stay + 90 Days)" , size(large)) xtitle("Differential Distance to PE NH (Miles)" , size(large)) line(qfit) note("Slope (per 10 Miles): `b' (`se')") n(10) ///
savegraph("$Output/Figures/4_Binscatter_Mort_Dist.png") replace

// BS 2
reg pay_res dist_res
local b = string(_b[dist_res]*10,"%9.3f")
local se = string(_se[dist_res]*10,"%9.3f")
binscatter pay_res dist_res, ytitle("Log Amount Billed per Patient Stay" , size(large)) xtitle("Differential Distance to PE NH (Miles)" , size(large)) line(qfit) note("Slope (per 10 Miles): `b' (`se')") n(10) ///
savegraph("$Output/Figures/4_Binscatter_Pay_Dist.png") replace

// BS 3
reg pay_res_post dist_res
local b = string(_b[dist_res]*10,"%9.3f")
local se = string(_se[dist_res]*10,"%9.3f")
binscatter pay_res_post dist_res, ytitle("Log Amount per Patient Stay + 90 Days" , size(large)) xtitle("Differential Distance to PE NH (Miles)" , size(large)) line(qfit) note("Slope (per 10 Miles): `b' (`se')") n(10) ///
savegraph("$Output/Figures/4_Binscatter_Pay_Post_Dist.png") replace

*************************************************
* MTE Analysis - Figure 5
*************************************************

*set parameters;
clear all
set maxvar 15000
*number of reps for bootstrap;
local reps 50
*mte polynomial order;
local order 2
*local min. bin count threshold 
local binthres 50
*90-day follow up period;
local durn 90
*90% confidence intervals;
local t =1.645


* Define instrument and covariate locals;

local gp gp_old_1 gp_old_2 gp_old_3 ///
	gp_old_4 gp_old_5 gp_old_6 gp_old_7 gp_old_8 gp_old_9 gp_old_10 gp_old_11 ///
	gp_old_12 gp_old_13 gp_old_14 gp_old_15 gp_old_16 gp_old_17

local z diff_dist_mod diff_dist_sq

local x female married dual black age female_age `gp' 

local varlist bene_id sex age charles_old dual married race black facid year peind logbeds ///
	bene_hrr diff_dist_mod diff_dist_sq mortality_30 mortality_90 gp_old* everpe ///
	HH_income_median_ip log_pmt length_stay_w log_length_stay


use `varlist' using "$Data/Temp_Files/Person_Base", clear
gen dd = diff_dist_mod*10

*Create controls;
gen female = sex==2
gen female_age = female*age

xtile dd_tile = dd, nq(100)

*Gen interactions;
local gp gp_old_1 gp_old_2 gp_old_3 ///
	gp_old_4 gp_old_5 gp_old_6 gp_old_7 gp_old_8 gp_old_9 gp_old_10 gp_old_11 ///
	gp_old_12 gp_old_13 gp_old_14 gp_old_15 gp_old_16 gp_old_17
	
local x female age female_age black dual married `gp'

foreach xvar of local x{		
	gen z_`xvar' = 	diff_dist_mod * `xvar'
	gen z2_`xvar' = diff_dist_sq * `xvar'
}

local z_gp z_gp_old_1 z_gp_old_2 z_gp_old_3 z_gp_old_4 z_gp_old_5 z_gp_old_6 z_gp_old_7 ///
	z_gp_old_8 z_gp_old_9 z_gp_old_10 z_gp_old_11 z_gp_old_12 z_gp_old_13 z_gp_old_14 ///
	z_gp_old_15 z_gp_old_16 z_gp_old_17
local z2_gp z2_gp_old_1 z2_gp_old_2 z2_gp_old_3 z2_gp_old_4 z2_gp_old_5 z2_gp_old_6 z2_gp_old_7 ///
	z2_gp_old_8 z2_gp_old_9 z2_gp_old_10 z2_gp_old_11 z2_gp_old_12 z2_gp_old_13 z2_gp_old_14 ///
	z2_gp_old_15 z2_gp_old_16 z2_gp_old_17	
local z_x z_female z_married z_dual z_black z_age z_female_age `z_gp'  ///
	  z2_female z2_married z2_dual z2_black z2_age z2_female_age `z2_gp' 
	  
local p_gp p_gp_old_1 p_gp_old_2 p_gp_old_3 ///
	p_gp_old_4 p_gp_old_5 p_gp_old_6 p_gp_old_7 p_gp_old_8 p_gp_old_9 p_gp_old_10 p_gp_old_11 ///
	p_gp_old_12 p_gp_old_13 p_gp_old_14 p_gp_old_15 p_gp_old_16 p_gp_old_17	

local p_x p_female p_married p_dual p_black p_age p_female_age `p_gp'
	  
***Run first stage; 	 
regress peind `x' `z_x' `z' i.facid bene_hrr#year, cluster(facid)
predict pscore

*Bound the predicted propensity to be in (0,1)
replace pscore = min(pscore,1)
replace pscore = max(pscore,0)

summ pscore
gen pnorm = r(mean)

gen pround = round(pscore,0.01)
replace pround = int(pround*100)/100


di "produce fig 5a - treatment propensity vs. differential distance"
preserve
	gen one =1
	collapse (sum) tot=one (mean) pscore dd, by(dd_tile)
	egen samp = total(tot)
	twoway (fpfit pscore dd, lc(navy) estopts(degree(4))) ///
	       (scatter pscore dd, msym(oh) mfc(white) mlc(navy)), ///	
		graphregion(color(white)) yla(0(0.1)0.6, format(%9.1f) angle(0)) legend(off) ///
		ytitle("P(z)") xla(-20(10)70) xtick(-20(5)70) xtitle("Differential distance (miles)")
	graph export "$Output/Figures/5a_Treatment_Diff_Distance.pdf", replace
restore

di "produce fig 5b - Common support figure"
preserve 
	gen one=1
	collapse (sum) tot_ = one , by(peind pround)
	bys peind: egen totcat_ = total(tot)
	reshape wide tot_ totcat_, i(pround) j(peind)
	keep if inrange(pround,0.01,0.99)

	twoway (bar tot_0 pround , lc(white) fcolor(red) color(%60) barwidth(0.01)) ///
	       (bar tot_1 pround , lc(navy) fc(white) yaxis(2) color(%60) barwidth(0.01)), ///
	       xtitle("Propensity Score") ytitle("Number of non-PE", axis(1)) ytitle("Number of PE", axis(2)) ///
	       legend(order(1 "Non-PE" 2 "PE")) xla(0(0.2)1) xtick(0(0.1)1) ///
	       graphregion(color(white)) yla(500 10000 100000 500000, angle(0) axis(1) format(%9.0gc) nogrid)  ///
	       yscale(log axis(1)) yla(50 1000 5000 40000, angle(0) axis(2) format(%8.0gc)) yscale(log axis(2)) 
	graph export "$Output/Figures/5b_Common_Support.pdf", replace     
	
restore 

* Define instrument and covariate locals;

local gp gp_old_1 gp_old_2 gp_old_3 ///
	gp_old_4 gp_old_5 gp_old_6 gp_old_7 gp_old_8 gp_old_9 gp_old_10 gp_old_11 ///
	gp_old_12 gp_old_13 gp_old_14 gp_old_15 gp_old_16 gp_old_17

local z diff_dist_mod diff_dist_sq

local x female married dual black age female_age `gp' 

local varlist bene_id sex age charles_old dual married race black facid year peind logbeds ///
	bene_hrr diff_dist_mod diff_dist_sq mortality_30 mortality_90 gp_old* everpe ///
	HH_income_median_ip log_pmt length_stay_w log_length_stay
	
*set parameters;
*number of reps for bootstrap;
local reps 50
*mte polynomial order;
local order 2
*local min. bin count threshold 
local binthres 50
*90-day follow up period;
local durn 90
*90% confidence intervals;
local t =1.645	

* Generate the IV effect and IV weights:
*Obtain residual of pscore after partialling out X;
regress pscore `x' i.facid bene_hrr#year
predict z_res if e(sample), res
tempfile mte
save `mte', replace

tempvar sample

*identify bins where we have fewer than 50 of each treated and control patients.
egen bincount = count(pscore), by(pround peind)
tab pround if bincount<`binthres' & peind==0
tab pround if bincount<`binthres' & peind==1

*Drop pscore bins p=(0.95,0.99) since they have fewer than 50 control patients;
drop if pround>0.94

*accordingly, set number of bins for unobserved component
local nbins 94

ivreghdfe mortality_`durn' `x' (peind = `z') , absorb(facid bene_hrr#year) tol(1e-6) accel(sd) cluster(facid)
local b_2sls  = round(_b[peind],0.001)
local se_2sls = _se[peind]
gen `sample' = e(sample)

sum peind if `sample'==1
local mean_d = r(mean)

sum z_res if `sample'==1
local mean_iv = r(mean)
local N=r(N)

qui gen cov = (z_res-`mean_iv')*(peind - `mean_d') if `sample'==1
sum cov
replace cov = cov/r(mean)

corr peind z_res if `sample'==1, cov
local cov = r(cov_12)
di `cov'

mat define weights_iv = J(1,`nbins',.)

foreach i of numlist 1/`nbins'{
	sum z_res if `sample'==1 & (pscore > (`i')/100)
	if r(N)>0 {
		mat weights_iv[1,`i'] = (r(mean) - `mean_iv')*r(N)/`N'/`nbins'/`cov'
	}	
	else{	
		mat weights_iv[1,`i'] = 0
	}
}

mat ivweight_sum = weights_iv * J(`nbins',1,1)
	
mat weights_iv = weights_iv / ivweight_sum[1,1]
mat list weights_iv

*Outcome equation using prop score;
*Polynomial in p;
 
local poly ""
forval p=2/`order'{
	gen p`p' = pscore^`p'
	local poly `poly' p`p'
}

di "polynomial is `poly'"

foreach var of varlist female married dual black age female_age gp_old* {
	gen p_`var' = pscore* `var'	
}

*Run outcome model;
tempfile mte_bs
bootstrap, reps(`reps') saving(`mte_bs', replace): reghdfe mortality_`durn' `x' `p_x' pscore `poly', absorb(facid bene_hrr#year) tol(1e-6) accel(sd) cluster(facid)

*Test of non-linearity of outcome w.r.t. propensity score i.e. treatment effects interact with selection;
test `poly'

*Bring in bootstrap estimates;
preserve 
	use `mte_bs', clear
	foreach var of varlist _all{	
		mkmat `var', matrix(m`var')
	}
restore

* Compute MTE and various weighted average TEs:	
*Compute the portion that is dependent on var values and coeffs. evaluated at mean Xs.
local mte_base = 0
matrix define mte_base_ci = J(`reps',1,0)

gen mte_x =0

foreach var of local p_x {	
	qui summ `var' if e(sample)
	local mte_base = `mte_base' + _b[`var'] * r(mean)
	replace mte_x = mte_x + _b[`var'] * `var'
}

*Confidence intervals for MTE, ATE, ATT, ATUT, and IV;
forval r=1/`reps'{
	gen mte_x`r' = 0	
	foreach var of local p_x {	
		qui summ `var' if e(sample)
		matrix mte_base_ci[`r',1] = mte_base_ci[`r',1] + m_b_`var'[`r',1] * r(mean)
		replace mte_x`r' = mte_x`r' + m_b_`var'[`r',1] * `var'
	}	
}

di "Observed part of MTE evaluated at means of Xs: " `mte_base'

qui summ mte_x 
local ate_base = r(mean)
di "Observed part of ATE evaluated at means of Xs: " `ate_base'

*Confidence intervals for ATE
matrix define ate_base_ci = J(`reps',1,.)
forval r = 1/`reps'{
	summ mte_x`r'
	matrix ate_base_ci[`r',1] = r(mean)
}

*Compute base value for TT
qui gen mtebase_tt = mte_x * pscore

*what does mtebase_tt look like by pscore
tab pround, summ(mtebase_tt)

*fix to be evaluated at X means
sum mtebase_tt
replace mtebase_tt = r(mean) if mtebase_tt<.
sum pscore
replace mtebase_tt = mtebase_tt / r(mean)
sum mtebase_tt
local tt_base = r(mean)

*Confidence interval for ATT
matrix define tt_base_ci = J(`reps',1,.)
forval r = 1/`reps'{
	gen mtebase_tt`r'= mte_x`r' * pscore
	summ mtebase_tt`r'
	replace mtebase_tt`r' = r(mean) if mtebase_tt`r'<.
	sum pscore
	replace mtebase_tt`r' = mtebase_tt`r' / r(mean)
	sum mtebase_tt`r'
	matrix tt_base_ci[`r',1] = r(mean)
	
	*delete since no longer required
	drop mtebase_tt`r'
}

*Now do the same for TUT;
gen mtebase_tut = mte_x * (1-pscore)

*what does mtebase_tut look like by pscore
tab pround, summ(mtebase_tut)

*fix to be evaluated at X means
sum mtebase_tut
replace mtebase_tut = r(mean) if mtebase_tut<.
sum pscore
replace mtebase_tut = mtebase_tut / (1-r(mean))
sum mtebase_tut
local tut_base = r(mean)

*Confidence intervals for ATUT;
matrix define tut_base_ci = J(`reps',1,.)
forval r = 1/`reps'{
	gen mtebase_tut`r'= mte_x`r' * (1-pscore)
	summ mtebase_tut`r'
	replace mtebase_tut`r' = r(mean) if mtebase_tut`r'<.
	sum pscore
	replace mtebase_tut`r' = mtebase_tut`r' / (1-r(mean))
	sum mtebase_tut`r'
	matrix tut_base_ci[`r',1] = r(mean)
	
	*delete since no longer required
	drop mtebase_tut`r'
}

*Compute base portion of IV at means of Xs;
gen mtebase_iv = mte_x * cov if pscore <.
sum mtebase_iv
replace mtebase_iv = r(mean) if mtebase_iv<.
local iv_base = r(mean)

*Confidence interval for IV
matrix define iv_base_ci = J(`reps',1,.)
forval r=1/`reps'{
	gen mtebase_iv`r' = mte_x`r' * cov if pscore<.
	sum mtebase_iv`r'
	replace mtebase_iv`r' = r(mean) if mtebase_iv`r' <.
	matrix iv_base_ci[`r',1] = r(mean)	
	drop mte_x`r'
}

*Combine base values with the control function by grid of u_d;

qui sum pscore
local mean_p = r(mean)
local totn_tt = 100*r(N)*r(mean)
local totn_tut = 100*r(N)*(1-r(mean))
	
matrix define mte = J(1,`nbins',.)
matrix define mte_iv = J(1,`nbins',.)
matrix define mte_unobs = J(1,`nbins',.)
matrix define tt_unobs_wt = J(1,`nbins',.)
matrix define tut_unobs_wt = J(1,`nbins',.)

forval i=1(1)`nbins'{	
	local mtePoly`i' = _b[pscore]
	forvalues j = 2(1)`order' {
			local mtePoly`i' = `mtePoly`i'' + `j' * _b[p`j'] * (`i' / 100) ^ (`j' - 1)
	}	
	 
	
	matrix mte[1,`i'] = `mte_base' + `mtePoly`i''
	matrix mte_unobs[1,`i'] = `mtePoly`i''
	
	matrix mte_iv[1,`i'] = `iv_base' + `mtePoly`i''
	
	qui sum pscore if pscore >`i'/100
	matrix tt_unobs_wt[1,`i'] = r(N)/`totn_tt'
	
	qui sum pscore if pscore <= `i'/100
	matrix tut_unobs_wt[1,`i'] = r(N)/`totn_tut'
}

*Confidence intervals;
matrix define mte_ci = J(`reps',`nbins',.)
matrix define mte_iv_ci = J(`reps',`nbins',.)
matrix define mte_unobs_ci = J(`reps',`nbins',.)

forval r=1/`reps'{
	forval i=1(1)`nbins'{	
		local mtePoly`i' = m_b_pscore[`r',1]
		forvalues j = 2(1)`order' {
			local mtePoly`i' = `mtePoly`i'' + `j' * m_b_p`j'[`r',1] * (`i' / 100) ^ (`j' - 1)
		}
	
		matrix mte_ci[`r',`i'] = mte_base_ci[`r',1] + `mtePoly`i''
		matrix mte_iv_ci[`r',`i'] = iv_base_ci[`r',1] + `mtePoly`i''
		matrix mte_unobs_ci[`r',`i'] = `mtePoly`i''
		
	}
}

*Compute ATE;
matrix define one = J(`nbins',1,1)
matrix ate_unobs = 1/`nbins' * (mte_unobs * one) 
matrix ate_unobs_ci  = 1/`nbins' * (mte_unobs_ci * one)
 
gen ate = `ate_base' + ate_unobs[1,1] in 1/`nbins'
summ ate 
local ate_mean = r(mean)
matrix ate_ci = ate_base_ci + ate_unobs_ci
svmat ate_ci
summ ate_ci1
local ate_lb = `ate_mean' - `t' * r(sd)
local ate_ub = `ate_mean' + `t' * r(sd)

di "ATE:  `ate_mean' (`ate_lb', `ate_ub')"
local ate_mean = round(`ate_mean',0.001)

*Compute TT and TUT;
matrix tt_unobs = mte_unobs * tt_unobs_wt'
gen tt = `tt_base' + tt_unobs[1,1] in 1/`nbins'
summ tt
local tt_mean = r(mean)

*now compute CI;
matrix tt_unobs_ci = mte_unobs_ci * tt_unobs_wt'
matrix tt_ci = tt_base_ci + tt_unobs_ci 
svmat tt_ci
summ tt_ci1
local tt_lb = `tt_mean' - `t' * r(sd)
local tt_ub = `tt_mean' + `t' * r(sd)

di "*************"
di "ATT: `tt_mean' (`tt_lb', `tt_ub')"
di "*************"
local tt_mean = round(`tt_mean',0.001)

matrix tut_unobs = mte_unobs * tut_unobs_wt'
gen tut = `tut_base' + tut_unobs[1,1] in 1/`nbins'
summ tut
local tut_mean = r(mean)

*now compute CI;
matrix tut_unobs_ci = mte_unobs_ci * tut_unobs_wt'
matrix tut_ci = tut_base_ci + tut_unobs_ci 
svmat tut_ci
summ tut_ci1
local tut_lb = `tut_mean' - `t' * r(sd)
local tut_ub = `tut_mean' + `t' * r(sd)

di "*************"
di "ATUT: `tut_mean' (`tut_lb' `tut_ub')"
di "*************"
local tut_mean = round(`tut_mean',0.001)

*Compute IV using MTE;
matrix iv = weights_iv * mte_iv'
gen iv = iv[1,1] in 1/`nbins'
summ iv
local iv_mean = r(mean)

matrix iv_ci = mte_iv_ci * weights_iv'
svmat iv_ci
summ iv_ci1
local iv_lb = `iv_mean' - `t' * r(sd)
local iv_ub = `iv_mean' + `t' * r(sd)

di "*************"
di "Weighted IV: `iv_mean' (`iv_lb' `iv_ub')"
di "*************"
local iv_mean = round(`iv_mean',0.001)

matrix b = mte'
svmat b
gen obs = _n in 1/`nbins'
gen ud = _n/100 in 1/`nbins'

*get ci of MTE values;
gen mte_sd=.


forval i=1/`nbins'{
	matrix mte`i' = mte_ci[1..`reps',`i']
	svmat mte`i'
	summ mte`i'1
	replace mte_sd = r(sd) in `i'
	matrix drop mte`i'
	drop mte`i'1
}

gen mte_lb = b1 - `t' * mte_sd
gen mte_ub = b1 + `t' * mte_sd

*Result figure:
*Fig: MTE curve with CI;
twoway (line b1 ud, lc(black)) ///
	(line mte_lb ud, lc(gs10) lp(dash)) ///
	(line ate ud, lc(red) lp(dash)) ///	
	(line mte_ub ud, lc(gs10) lp(dash)), ///
	graphregion(color(white)) ///
	xtitle("Unobserved resistance to treatment") xla(0(0.1)1, format(%9.1f)) xtick(0(0.05)1)  ///
	ytitle("MTE") yla(-0.04(0.01)0.08, angle(0) nogrid format(%9.2f))  ///
	legend(order(1 "MTE" 2 "90% CI" 3 "ATE (`ate_mean')") rows(1))
graph export "$Output/Figures/5c_MTE_CI.pdf", replace

*Fig 2: TT and TUT weights with ATT and ATUT lines;

mat ttwt = tt_unobs_wt'
svmat ttwt 
mat tutwt = tut_unobs_wt'
svmat tutwt 

twoway (scatter ttwt1 ud, msymb(O) mfc(gs5) mlc(gs5) yaxis(2)) ///
	(scatter tutwt1 ud, msymb(T) mfc(gs10) mlc(gs10) yaxis(2)) ///
	(line tt ud, lc(gs5) lp(longdash) yaxis(1)) ///
	(line tut ud, lc(gs10) lp(shortdash) yaxis(1)), ///
	graphregion(color(white)) ///
	xtitle("Unobserved resistance to treatment") xla(0(0.1)1, format(%9.1f)) xtick(0(0.05)1)  ///
	ytitle("Treatment effect", axis(1)) yla(0(0.01)0.05, angle(0) nogrid axis(1) format(%9.3f))  ///
	ytitle("Weights", axis(2)) yla(0(0.01)0.05, angle(0) nogrid axis(2) format(%9.3f))  ///
	legend(order(1 "TT weights" 2 "TUT weights" 3 "ATT (`tt_mean')" 4 "ATUT (`tut_mean')") rows(2))
graph export "$Output/Figures/5d_TT.pdf", replace


*************************************************
* Facility Event Studies - Figures 6 and 7
*************************************************
use "$Data/Temp_Files/Facility_Base", replace
// Generate dummies
destring deal_year, replace
gen event_time=year-deal_year
replace event_time = 0 if mi(event_time)
tab event_time, gen(deal_dummies)

label var rating_overall_w "Overall Rating"
label var rating_survey_w  "Deficiency Rating"
label var staff_hr_by_resday_w "Staff Hours Per Resident Day"
label var cna_hr_by_resday_w "Nurse Assistant Hours Per Resident Day"
label var lpn_hr_by_resday_w "Licensed Nurse Hours Per Resident Day"
label var rn_hr_by_resday_w "Registered nurse Hours Per Resident Day"
label var lman_fee_all_totalxw "Log Management Fee Cost"
label var llease_all_totalxw "Log Building Lease Cost"
label var lint_expense_all_totalxw "Log Interest Cost"


local dummies deal_dummies1-deal_dummies4 deal_dummies6-deal_dummies11
local outcomes rating_overall_w rating_survey_w  staff_hr_by_resday_w cna_hr_by_resday_w lpn_hr_by_resday_w rn_hr_by_resday_w  
// Figure 6
foreach depvar of local outcomes {
	local vtext: variable label `depvar'
	di "`depvar'"
	* Regression
	reghdfe `depvar' `dummies' la_herfbeds_cty res_profit_cty res_multifac_cty res_hospbase_cty logcount agg_cmi_w acuityindx_w agg_black_w , absorb(year facid) cluster(facid)

	* Save for figure
	local var outcome
	gen `var'_b = .
	gen `var'_sd = .
	gen `var'_u = .
	gen `var'_l = .
	
	* pick years before; with 17 as excluded  
	forvalues x=2(1)4{
		local y = `x' - 1
		replace `var'_b=_b[deal_dummies`x'] in `y'
		replace `var'_sd=_se[deal_dummies`x'] in `y'
		replace `var'_u=_b[deal_dummies`x']+1.96*_se[deal_dummies`x'] in `y'
		replace `var'_l=_b[deal_dummies`x']-1.96*_se[deal_dummies`x'] in `y'
	}
	
	* Baseline
	replace `var'_b = 0 in 4
	
	* fill in for years after 
	* UPDATE if data changes (i.e. deal year == 0 dummy can change)
	forvalues x=6(1)10{
		local y = `x' - 1
		replace `var'_b=_b[deal_dummies`x'] in `y'
		replace `var'_sd=_se[deal_dummies`x'] in `y'
		replace `var'_u=_b[deal_dummies`x']+1.96*_se[deal_dummies`x'] in `y'
		replace `var'_l=_b[deal_dummies`x']-1.96*_se[deal_dummies`x'] in `y'
	}
	
	* Figure label cleanup
	gen n=_n
	label define betas   1 "-4" 2 "-3" 3 "-2" 4 "-1" 5 "0" 6 "1" 7 "2" 8 "3"  9 "4" , replace 
	label values n betas
		
	
	twoway (rcap `var'_u `var'_l n, color(gs11)) (scatter `var'_b n, msize(medium) mcolor(navy) msymbol(circle) )  if n <= 9,  graphregion(color(white)) legend(off)  xlabel(1(1)9,valuelabel) ytitle("`vtext'", size(medsmall))  xtitle("Years from Deal (0 is Deal Year)")   yline(0, lcolor(black) ) xline(4.5, lcolor(red) lpattern(dash)) 
	graph export "$Output/Figures/6_`depvar'.png", replace 
 
		
	drop `var'_b
	drop `var'_sd
	drop `var'_u
	drop `var'_l
	drop n 
}

local outcomes lman_fee_all_totalxw llease_all_totalxw lint_expense_all_totalxw
// Figure 7
foreach depvar of local outcomes {
	local vtext: variable label `depvar'
	di "`depvar'"
	* Regression
	reghdfe `depvar' `dummies' la_herfbeds_cty res_profit_cty res_multifac_cty res_hospbase_cty logcount agg_cmi_w acuityindx_w agg_black_w , absorb(year facid) cluster(facid)

	* Save for figure
	local var outcome
	gen `var'_b = .
	gen `var'_sd = .
	gen `var'_u = .
	gen `var'_l = .
	
	* pick years before; with 17 as excluded  
	forvalues x=2(1)4{
		local y = `x' - 1
		replace `var'_b=_b[deal_dummies`x'] in `y'
		replace `var'_sd=_se[deal_dummies`x'] in `y'
		replace `var'_u=_b[deal_dummies`x']+1.96*_se[deal_dummies`x'] in `y'
		replace `var'_l=_b[deal_dummies`x']-1.96*_se[deal_dummies`x'] in `y'
	}
	
	* Baseline
	replace `var'_b = 0 in 4
	
	* fill in for years after 
	* UPDATE if data changes (i.e. deal year == 0 dummy can change)
	forvalues x=6(1)10{
		local y = `x' - 1
		replace `var'_b=_b[deal_dummies`x'] in `y'
		replace `var'_sd=_se[deal_dummies`x'] in `y'
		replace `var'_u=_b[deal_dummies`x']+1.96*_se[deal_dummies`x'] in `y'
		replace `var'_l=_b[deal_dummies`x']-1.96*_se[deal_dummies`x'] in `y'
	}
	
	* Figure label cleanup
	gen n=_n
	label define betas   1 "-4" 2 "-3" 3 "-2" 4 "-1" 5 "0" 6 "1" 7 "2" 8 "3"  9 "4" , replace 
	label values n betas
		
	
	twoway (rcap `var'_u `var'_l n, color(gs11)) (scatter `var'_b n, msize(medium) mcolor(navy) msymbol(circle) )  if n <= 9,  graphregion(color(white)) legend(off)  xlabel(1(1)9,valuelabel) ytitle("`vtext'", size(medsmall))  xtitle("Years from Deal (0 is Deal Year)")   yline(0, lcolor(black) ) xline(4.5, lcolor(red) lpattern(dash)) 
	graph export "$Output/Figures/7_`depvar'.png", replace 
 
		
	drop `var'_b
	drop `var'_sd
	drop `var'_u
	drop `var'_l
	drop n 
}